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ABSTRACT 



CO 

o 

■ The mechanisms which could lead to chemo-thermal instabilities and fragmentation 

during the formation of primordial protostars are investigated analytically. We in- 
troduce new analytic approximations for H2 cooling rates bridging the optically thin 
and thick regimes. These allow us to discuss chemo-thermal instabilities up to densities 
when protostars become optically thick to continuum radiation (n = p/mn ^ 10 16-17 ci 

I/"*) I During the proto-stellar collapse instabilities are active in two different density regimes. 

In the well known "low density" regime (n ~ 10 s — 10 10 cm -3 ), instability is due to 
3-body reactions quickly converting atomic hydrogen into H 2 . In the "high density" 

■ regime (n £ 10 cm -3 ), another instability is triggered by the strong increase in the 
' cooling rate due to H2 Collisional Induced Emission (CIE). In agreement with the 

■ three dimensional simulations, we find that the "low density" instabilities cannot lead 

to fragmentation, both because fluctuations are too small to survive turbulent mix- 
ing, and because their growth times are too slow. The situation for the newly found 

I "high density" instability is analytically similar. This continuum cooling instability 

is as weak as "low density" instability, with similar ratios of growth and dynamical 
time scales, as well as allowing for the necessary fragmentation condition t coo i ^ td yn . 
O Because the instability growth timescale is always longer than the free fall timescale, 

it seems unlikely that fragmentation could occur in this high density regime. Con- 
sequently, one expects the first stars to be very massive, not to form binaries nor 
^ ' harbour planets. Nevertheless, full three dimensional simulations are required to be 

certain. Such 3D calculations could become possible using simplified approaches to 
approximate the effects of radiative transfer, which we show to work very well in ID 
calculations, giving virtually indistinguishable results from calculations employing full 
y\ ' line transfer. This indicates that the effects of radiative transfer during the initial 

^ . stages of formation of primordial proto-stars are local corrections to cooling rather 

than influencing the energetics of distant regions of the flow. 

Key words: stars: formation - instabilities - molecular processes. 



1 INTRODUCTION 

In current scenarios for the formation of the first stars, molecular hydrogen plays a prominent role, providing the only 
cooling mechanism for metal-free gas at T ^ 10 4 K. H2 cooling is responsible for both the collapse of the gas inside the 
small (Afhaio ~ 10 6 M@) cosmological halos where primordial star formation is believed to occur {e.g. Couchman & Rees 
1986; Tegmark et al. 1997, Abel et al. 1998), and for the fragmentation of the gas itself into "clumps" with masses of 
~ 100 — 1000 M©, as is seen in the simulations of Abel, Bryan & Norman (1998, 2000) (see also Bromm, Coppi & Larson 
2002). The fate of these objects is unclear. Several previous studies, based on analytical arguments (i.e. stability analysis) 
or single-zone models, led to different conclusions about the fragmentation properties of primordial clouds. One important 
instability-triggering mechanism was first suggested by Palla, Salpeter & Stahler (1983; hereafter PSS83): the onset of 3-body 
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H2 formation at number densities n Si 10 8 cm -3 (where n = p/mn) causes a fast increase in the cooling rate, possibly leading 
to fragmentation on a mass scale of ~ O.IA/q. A similar result was obtained by Silk (1983; hereafter S83) by applying a 
criterion for "chemo-thermal instability" (first discussed by Sabano & Yoshii, 1977) to an object where 3-body H2 formation 
is important. 

Numerical simulations have recently shed new light upon the issue, and in particular Abel, Bryan & Norman (2002; 
hereafter ABN02) were able to reach the 3-body reactions density regime finding that these reactions actually lead to the 
formation of a single fully molecular core (with a mass of ~ 1 Mq) at the centre of each clump. They point out that in their 
simulations they see several chemo-thermally unstable regions, but that the instabilities do not lead to fragmentation because 
turbulence efficiently mixes the gas, erasing each fluctuation before it can grow significantly. Recently, Omukai & Yoshii (2003; 
OY03) found a somewhat different explanation by improving the original S83 investigation and pointing out that instabilities 
growth was too slow. 

The further evolution of such an object is still uncertain: due to the high molecular fraction, many of the numerous 
H2 roto-vibrational lines (accounting for most of the cooling) become optically thick just after the formation of the molecular 
core. For this reason, it becomes impossible to estimate the H2 cooling rate without a complete treatment of radiative transfer, 
which currently can not be included into full 3-D hydrodynamical simulations. Simpler 1-D studies, such as those in Omukai 
& Nishi (1998; hereafter ON98) and Ripamonti et al. (2002; hereafter R02), can follow the further evolution, giving useful 
predictions about the physical processes driving the collapse, or about the properties of the small hydrostatic protostellar 
core which finally forms in the centre; on the other hand, their intrinsic spherical symmetry prevents them from giving direct 
predictions about fragmentation. 

An interesting result of ON98 and R02 is that during the collapse, the protostellar object experiences a phase when 
cooling is dominated by the effects of H2 Collision-Induced Emission (CIE), the opposite process of the more commonly known 
Collision-Induced Absorption, or CIA. Both of them find that, at the beginning of the CIE-cooling phase, the molecular core 
is optically thin to CIE-produced continuum photons, and there is a very fast increase in the cooling rate until the continuum 
optical depth exceeds unity. However, neither of these previous studies pointed out that during this phase the core conditions 
resemble the ones that, at the onset of 3-body H2 formation, led to the chemo-thermal instability: the core is optically thin 
(which is believed to be a necessary condition for fragmentation; see Rees 1976) and the cooling rate is undergoing a dramatic 
increase. 

In the next sections, we will investigate more thoroughly the issue of chemo-thermal instability, with special attention to 
the effects of CIE cooling. In section 2 we will describe the main cooling processes, giving a brief account of CIE emission, and 
some useful approximations for H2 line cooling. In section 3 we examine the conditions for the insurgence of chemo-thermal 
instabilities, and argue whether they can cause fragmentation. Finally, we discuss the results in section 4. 



2 COOLING PROCESSES 
2.1 Collision-Induced Emission 

H2 molecules have no electric dipole, and emission or absorption of radiation can take place only through quadrupole transi- 
tions. But when a collision takes place, the interacting pair (H2-H2, tb-He, H2-H) briefly acts as a "supermolecule" with a 
nonzero electric dipole, and an high probability of emitting (CIE) or absorbing (CIA) a photon (see the brief discussion in 
Lenzuni, Chernoff & Salpeter 1991, or Frommhold 1993 for an extensive account). In the same way, during an H-He collision, 
the two atoms perturb each other and there is a relevant probability of emission (or absorption) of a photon through a dipole 
transition. Because of the very short collision times (At 10~ 12 s for T Si 300 K), collision-induced lines become very broad, 
actually merging into a continuum; for example, in the H2 CIE spectrum only the vibrational bands can be discerned as 
smooth peaks deriving from the merging of the roto-vibrational lines. 



2.1.1 CIE cooling rate 

The shape of CIE spectra can be found in the literature (see Table Q for a list of references), and used for estimating the 
cooling rates due to the main CIE processes. The definition of the monochromatic emission coefficient j v (cfr. Rybicki & 
Lightman, eq. 1.15) 

dE = j v dV dQdtdu (1) 

(where E is the energy emitted, and j v has units of erg cm -3 ster -1 s _1 Hz -1 ) can be easily integrated to give the luminosity 
per unit mass L 

dE dE 47T f . , 

— / J.du. (2) 



dt dm dt pdV 
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Table 1. List of references for the various kind of collisions leading to CIE emission. We also specify the temperature range and the 
maximum frequency considered. 



Pair 


T[K] 


f max [em 1 ] 


Reference 


H 2 -H 2 


400-1000 


17000 


Borysow 2002 


H 2 -H 2 


1000-7000 


20000 


Borysow et al. 2001 


H 2 -He 


1000-7000 


20088 


Jorgensen et al. 2000 


H 2 -H 


1000-2500 


10000 


Gustafsson, Frommhold 2003 


H 2 -H 


400-1000 


6000 


Gustafsson et al. 2003 


H-He 


1500-10000 


11000 


Gustafsson, Frommhold 2001 



The gas can be assumed to be in thermal equilibrium, so that j v — a v B v (T) (where B V (T) is the Planck function at 
temperature T), and the cooling rate per unit mass due to emission induced by collisions between particles of species i and 
species j in a gas of temperature T and density p is 

47T f 

Lt,j{T,P,Vt,Vi) = — a„,i tj (T,P,Vi,Vi)B*{T)dv, (3) 

where yi and yj are the mass fractions of particles of species i and j, respectively. We took the values of the collision-induced 
absorption coefficient a Vl ij from the references listed in Table 0\ 

It can be seen that these references provide a complete set of data only in a relatively narrow temperature range 
(1500 K < T < 2500 K). By a fortunate coincidence, this is also the temperature range where CIE cooling dominates over all 
the other cooling mechanisms (see next section), and where CIE effects are more pronounced. However, we have artificially 
extended this range to the full extension of the H 2 -H 2 data range (400 K < T < 7000 K) by estimating the cooling rates 
Lij(T) which are not directly available, through the assumption that the ratio Lij(T) / 'Lh 2i h 2 (T) is the same as at the 
closest known value. This assumption is somewhat arbitrary, but has negligible effects for a gas where H 2 is the dominant 
species (which is the most relevant case); in the case of a gas with a high atomic fraction, and outside the 1500 K — 2500 K 
temperature range, the results shown in Fig. are only indicative. 

In Fig. Q we compare the CIE cooling rates of a pure H 2 gas at different densities (similar results hold also for H 2 -He and 
H 2 -H-He mixtures) with the H 2 line cooling rates calculated by several other authors. CIE cooling overcomes roto-vibrational 
line cooling at an H 2 number density tih 2 between 10 14 and 10 16 cm -3 , the exact value depending on the gas temperature. 
This values can change because of the effects of optical depth, and are appropriate only when both CIE cooling and H 2 line 
cooling occur in an optically thin regime. Instead, both ON98 and R02 (treating CIE by means of the Planck opacities given 
by Lenzuni et al. 1991) found that H 2 line cooling starts to be limited by optical depth at an early stage, when CIE cooling is 
still optically thin. As a result, in their models CIE starts to overcome H 2 lines at a relatively low density (corresponding to 
T c ~ 1600 K, riH 2 ,c — 5 x 10 13 cm -3 , where the c subscript refers to the conditions in their central regions; see Fig. 4 of R02). 

The gas chemical composition determines which kind of collisions contribute most to total cooling. In Fig. [5] we show the 
temperature dependency of the cooling rate due to each kind of collisions for three different molecular fractions, /h 2 = 0.999, 
0.5 and 0.001, where /h 2 = nH + "n H ' s ^ ne mass fraction of H 2 (relative to all the H atoms). The He mass fraction is kept fixed 
at j/ho = 0.25 (the value we will use throughout the rest of this paper), and the assumed density is p = 1.67 x 10 _10 gcm -3 
(i.e. n — 10 14 cm -3 ). As can be expected, H 2 collisions (especially H 2 -H 2 ) dominate CIE cooling when hydrogen is mostly 
molecular, but their importance declines with decreasing H 2 fraction, becoming negligibly small for mostly atomic gas. The 
total CIE cooling rate for molecular gas is substantially (20-30 times) larger than for atomic gas in the same conditions. 

2.1.2 Approximate CIE cooling rate 

In the rest of this paper, the only relevant case will be the one with /h 2 — 1. In that case, the total CIE cooling rate (per 
unit mass) 

Lcte{p, T, j/h, yn 2 , 3/Hc) = £h 2 ,h 2 + £h 2 ,h + £h 2 ,h g + Ln,iie (4) 
can be reasonably approximated by a power law 

Lcmip,T,yH,yn 2 ,yne) ~ Lcie (p, F, /h 2 ) = A C ie pT a Xfn 2 (5) 



1 These references actually give the absorption coefficients in cm -1 amagat -2 , which must be multiplied by a factor (ptji/ po,i)(pyj /poj) 
in order to find the absorption coefficients appropriate for the given density and mass fractions. The densities po,i are given by po,i = 
(Po?t),;)/(A;bTo), where Pq = 1 atm, To = 273.15 K, mi is the mass of a particle of species i and fee i s the Boltzmann constant. 
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Figure 1. Comparison of H2 line cooling rate and H2-H2 CIE cooling rate. Thick lines denote the line cooling rate as calculated by Galli 
& Palla 1998 (short dashes), Martin, Schwarz & Mandy 1996 (long dashes) and Lepp & Shull 1983 (dots); these were calculated assuming 
to be in the high density regime (n zi 10 8 cm -3 ), where the cooling per molecule is independent of density. Thin solid lines denote CIE 
cooling rates at molecular number densities nn 2 = 10 14 cm -3 (bottom line), n,fj 2 = 10 15 cm -3 (middle line) and nn 2 = 10 16 cm -3 (top 
line). If the considered object were optically thin to H2 line cooling, CIE would become important at number densities n ^ 10 15 cm -3 . 

with a ~ 4.0 and Acie — 0.072 erg cm 3 g -2 s _1 K _4 '°, and where X = 1 — j/h c = 0.75 is the total hydrogen mass fraction. As 
can be seen in fig.|3] this approximation holds down to /h 2 ~ 0.5, although very roughly. 



2.2 H2 line cooling 

2.2.1 Optically thin H2 line cooling 

At the moderately high densities we are interested in (n ^ 10 6 cm -3 ), the populations of H2 roto- vibrational levels can be 
found through the LTE assumption, and the H2 line cooling is quite well known (see Fig. 0. Here we will use the LTE rate 
first given by Hollenbach & McKee (1979), which was also used by Galli & Palla (1998) as an high density limit. 



-Klines, HM(p, T, /h 2 ) 



ran 



H(T 3 



(6) 



where Limes, hm is the H2 lines cooling rate per unit mass, T3 

( 0-13 'V' 3 

H(T 3 )-- i ^ — i e v~) 



and 



1+0. 12T^ 



6.7 x 10" 19 e T 3 + 1.6 x 10~ 18 



0.51 
" T 3 + 



11.7 

e T 3 



ergs 



(7) 



Eq. © is actually an analytic fit to the sum of the luminosities of all the roto-vibrational lines (which instead was used 
by ON98 and R02): 



Llixs 



2mu ^ 



hVm.nA 



2J n + 1 



U(T) 



(8) 



where v m ,n is the frequency of the H2 transition from level n downward to level m, A m ,n is the spontaneous radiative decay 
rate (Einstein coefficient) of the transition, J n is the rotational quantum number associated with level n, E n is the energy of 
level n and U(T) = (2Ji + l)e Bi// ' fcBT ' ) is the H2 partition function at temperature T. The sum in JHJ is over all H2 lines, 
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Figure 2. Comparison of CIE cooling rates for different kinds of collisions. Top panels show the cooling rates per unit mass as a function 
of temperature, for a gas with different molecular fractions (/h 2 = 0.999, 0.5 and 0.001). The thick solid line denotes the total CIE 
cooling, the thin lines denote the single components: H2-H2 (thin solid), H2-He (short dashed), H2-H (dot dashed) and H-He (long 
dashed); the dotted lines in the upper left and upper center panels show the results of eq. JFJ- Note that in the top-left panel the total 
and H2-H2 lines are practically indistinguishable, as well as the total and the H-He lines in the top-right panel. The bottom panels 
show the ratios of the last three components to the CIE cooling rate due to H2-H2 collisions; in these panels the line thickness indicates 
whether the ratio was taken from one of the references in table Q (thick lines) or it was assumed to be constant at the closest known 
value (thin lines). All the quantities were calculated assuming n = 10 14 cm -3 and j/hc = 0.25. 

with each term representing the total luminosity of an H2 molecule in the considered line. We note that the term in square 
brackets represents the fraction of H2 molecules in the roto- vibrational level n. 



2.2.2 Optically thick H2 line cooling 

The above cooling rates (eq. El and can only be applied to optically thin objects. In reality, as 3-body reactions start 
transforming the bulk of hydrogen into molecular form, the core regions of the contracting protostellar clouds quickly become 
optically thick to H2 line radiation (ON98, R02). As a result, the effective cooling rate falls well below the one predicted in 
the optically thin case. 

It is possible to find a reasonable approximation to the "effective" cooling rate in the central regions of a collapsing 
protostellar object. The density profile of such an object can be approximated with a central flat "core" (where also temperature 
and chemical composition are approximately constant) surrounded by an "envelope" where the density decreases as a power- 
law (p oc r~ 2,2 , according to both ON98 and R02). Such a density profile is predicted by the Larson-Penston self-similar 
solution (Larson 1969; Penston 1969), which applies fairly well. Another correct prediction that can be inferred from the 
Larson-Penston solution is that the mass of the central "core" is of the order of the Bonnor-Ebert mass {i.e. the critical 
mass for gravitational collapse of an isothermal sphere given an external pressure P ox t; see Ebert 1955 and Bonnor 1956), 
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Table 2. Comparison of the core properties as found at several stages of R02 calculations, to the values employed in the text: mass 
and radius are compared to the Bonnor Ebert values (ea. l9l and 1101 . while the column densities (Nn 2iOC and Nn 2 , COIC , representing the 
column density from the centre to infinity and from the centre to the edge of the core, respectively) are normalized to c n c X fn 2 C /2. 
In order to extract core data from the simulations, we defined it as the central region where the density dependence on radius is flatter 
than p tx r~ ■ 



T C [K] 


n c [a.m.u. cm ] 


/h 2 ,c 


M core /AfBE,c 


-Rcore/ RbE.c 


JVH2,oo _ f 
H BE i c Tl c A JH 2 ,c/' i 


JVH 2 ,corc 
^BE,c n c^/H 2 ,c/2 


300 


9.1 X 10 6 


0.001 


0.24 


0.69 


0.62 


0.41 


401 


7.2 x 10 7 


0.002 


0.22 


0.67 


0.54 


0.41 


502 


1.0 x 10 9 


0.02 


0.20 


0.64 


0.45 


0.39 


594 


5.5 x 10 9 


0.11 


0.22 


0.67 


0.48 


0.42 


740 


2.8 x 10 10 


0.39 


0.29 


0.76 


0.60 


0.51 


913 


1.4 x 10 11 


0.72 


0.78 


1.17 


0.81 


0.62 


1328 


5.5 x 10 12 


0.98 


0.49 


0.88 


1.23 


0.72 


1636 


4.8 x 10 13 


0.98 


0.57 


0.95 


1.31 


0.73 


1799 


1.8 x 10 14 


0.96 


0.58 


0.94 


1.40 


0.79 


2007 


1.1 x 10 15 


0.92 


0.64 


0.99 


1.37 


0.74 


2200 


8.0 x 10 15 


0.91 


0.69 


1.00 


1.43 


0.77 


2337 


4.6 x 10 16 


0.92 


0.85 


1.12 


1.41 


0.77 



as calculated using the central value of temperature (T c ), number density (n c = p c /mn), molecular weight (fi c ) and mean 
adiabatic index (7 C ) 

/ rp x 3/2 , v -1/2 

A/be. c ^2OM ( t ^) J »ZH\ (9) 
from this a "Bonnor-Ebert radius" can be estimated 

„ /SMbe^V 73 1SJv1n i9 fT c \ l/2 f nc \~ 1/2 -2/3 2/3 , 1n * 

Rbe - = ita) * x 10 cm ItkJ It^J ^ 7c • (10) 

The H2 column density from the centre to infinity is then 

M ~ r ncXfH 2 ,t f°° n c Xf H2 {r) ( r \ 22 

N ** e - ,c 2 -/rbe.c — 2 — (jwj dr - 

n c A/H 2 ,c I . . / r , , „ w , 1-2.2, 



= Rbe,c * 2 ' |1 + J [fH 2 (xR3E, c )/fn 2 ,c]x dx^ (11) 
where x = t/Rbe,c- This last equation leads us to define a parameter £ such that 

iV H2 ,c = (flBE,c " J 2 /H2 ' C (12) 

and whose approximate value (cfr. eq. 111! is 

/oc 
[/h 2 (lflBE,c)//H 2|C ]l^ 2 ' 2 dl. (13) 

Since we are considering objects in which /h 2 is maximum at the centre, equation 11131 implies that 1 < £ < 1.8, but in 
Table|21(see columns 6 and 7) we show that this is not the case. The main reasons are that the density inside the core actually 
decreases from the centre outward, and that the Bonnor-Ebert mass provides an estimate of the core size which is not very 
accurate (Table |5J columns 4 and 5). In the following we will just set £ to values in agreement with the results of R02 (for 
example, Fig. [3]was obtained using £ = 0.20). 

Inside the "core", where we can assume that the temperature is reasonably constant, the cross section due to a transition 
from level m to level n is (see Lang 1980, eq. 2-69) 

c 2 / H^n \ 

<Jm, n {i>) = (f>(T c ,v,v m>n )- — = — e fc B^o - 1 A m> „, (14) 
07rfm,n \ J 

where <j> is the line profile. This cross section applies only to H2 molecules in the m roto-vibrational level, that is to a fraction 

Erg 

2 u{j^) e kBT " °f au t ne H2 molecules. 

In the range of conditions we are considering, the line profile is determined by Doppler broadening, with typical width 



Au D (T c ) = ^J^^. (15) 
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Figure 3. Ratio of cooling rates to the optically thin H2 line cooling rate, as a function of density. The solid thick lines show the evolution 
of the effective (that is, as influenced by optical depth) H2 line cooling rate as calculated by R02 in four different shells (the lowest one 
being the central shell), and the thick dot-dashed line shows the CIE cooling rate. The thin dashed, long dashed and dotted lines show 
the predictions of eq. 1191 when assuming T c = 500 K, 1000 K and 2000 K, respectively (note that these lines are not in monotonic order; 
the T c = 500 K line actually lies between the other two); instead, the thick dashed line shows the prediction of the same equation when 
the temperature evolution of the central shell is kept into account. Finally, the thin solid line shows the simple power law approximation 
of eq. |22j. All the curves making use ofeq. G2J assume £ = 0.20. 



Here, we will assume that the line profile can be simplified to 

^(T v v ) = { 2^d(t c ) if lf-fm.nl < Ai/d(T c ) ( "inside" the line) ; 

} if \v — v m ,n \ < Ai/d (Tc) ( "outside" the line), 



(16) 



and we will neglect the Doppler shifts due to the bulk motions 2 . With this approximation, and averaging over all H2 molecules, 
the mean cross section of a generic H2 molecule to a photon with frequency within Avu(T c ) from v m>n is 

„3 



m H 



1 Ar, 



2J m + 1 



k B T c \ J ' U(T C ) 

and the resulting optical depth from centre to infinity is then 

2J m + l e 



0.53 £ 



- 1 



U(T C 



A-T, 



10- 



2/3 



(17) 



(18) 



From this result, we can obtain an average cooling rate (per unit mass) inside the "core" by estimating the total luminosity 
of the core as seen from outside its surface. This can be done by assuming that the "core photosphere" (the region of the core 
where the optical depth to the surface is < 1) radiates in the optically thin limit, while the interior (where the optical depth 
to the surface is > 1) does not contribute to the total luminosity. 

Since the optical depths of the various lines are different, this amounts introducing correcting factors the for cooling rate 
of each line into eq. JHJ: 



hick(r) 



^7h 2 



hfrn.nA 



m.n-^m.n 



2J n + l 



U(T) 



(19) 



2 At the considered stages, the collapsing cloud is not in hydrostatic equilibrium, not even in the core: bulk velocities can amount to a 
few xlO 5 cms -1 . The resulting frequency shift is smaller than the line width, although not completely negligible. 
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Figure 4. Comparison of the results of a modified version of the R02 code (solid lines) and of the original one (dashed lines). In the 
modified version the cooling rate is given by eq. 1361 , which is the sum of the approximate cooling rates discussed in Section 2.2.2 
(H2 lines) and in Sections 2.2.1 and 3.3 (CIE); instead, in the original code the cooling rate is estimated through radiative transfer 
calculations. In the top panel we compare the evolution of four shells (including the central one) in the n — T plane; since all the lines 
were very close to each other, we artificially shifted down the results of the three non-central shells (with the bottom lines corresponding 
to the outermost shell). In the bottom panel we compare the evolution of the same four shells in the n — v plane. Note that two of our 
zones do not reach the maximum density shown in these plots, so the corresponding lines stop at n ~ 10 14 cm -3 and n ~ 10 16 cm -3 



where V m ,n is the ratio between the volume of the "core photosphere" for the considered (m, n) transition and the total core 
volume, given by 









1 - 


('- 


-i)1 





if T m ,n < 1 
if Tm.n > 1. 



while the geometrical correction 
f 1 



Grn,n 



1/2 



if T m , n < 1 
if Tm.n > 1. 



(20) 



(21) 



accounts for the radiation emitted towards the interior. 

We note that both the form of V m ,n as the ratio of the photospheric to the core volume, and the photospheric volume 
itself (at least in the optically thick case) clearly depend on the assumption that inside the core density and temperature are 
constant. 

In Fig.dlwe compare the results of eq. 1191 to the findings of R02. Provided that a sensible value of £ is chosen (in Fig. 01 
we are using £ = 0.20), the agreement is very good up to number densities of n ~ 10 13 cm -3 , and remains reasonable (within 
a factor of 2) up to the density where CIE cooling starts to be dominant (n ~ 10 14 cm -3 ). The differences can be partially 
explained by noting that at the highest densities the £ values reported in Table |2| are substantially higher than the adopted 
value of 0.20. It is also remarkable that the agreement remains good also when we consider the R02 results for non-central 
shells, provided that their own physical properties (T, n and /h 2 ) are plugged into eq. 1181 instead of the central ones. 
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Figure 5. Comparison of the results of the original and of the modified R02 code (see caption of fig. |1J. The various panels show the 
spatial profiles of density (top left), velocity (top right), temperature (bottom left) and molecular fraction (bottom right) at five different 
stages (specified through their central temperature) when the simplification in the H2 lines cooling rate is relevant. Dashed lines show 
the results of the original code, while solid lines show the result of the simplified one. The spatial coordinate is specified through the 
Lagrangian coordinate of R02, i.e. the enclosed mass. 



In the next sections we need an analytic expression for the cooling rate which is less cumbersome than eq. 1191 . We choose 
to employ a very simple approximation of the form 

Limes,thick(T) = Li incs , thin (T) min (1, (n/n ) _/3 ) (22) 

with n = 8 x 10 9 cm -3 (equivalent to p = 1.34 x 1(T 14 gem -3 ) and (3 = 0.45. 

We have tested the validity of both equations 1191 and 1221 by substituting it to the radiative transfer treatment inside 
the code used by R02 for their simulations, and comparing the results with those of the original code; in these tests we also 
added an approximate CIE contribution as can be found in eq. 1361 . As can be seen in Fig. 0] and |S] (which show the results 
obtained using eg. 1221 . the evolution is strikingly similar, not only in the central zones but also in the outer regions. So, we 
have been able to predict the global effect of radiative transfer on the H2 line cooling with a simple model based on purely local 
properties. It is quite remarkable (and partly unexpected) that a treatment which was believed to apply just to the centre of 
the "core" can be applied also to the non-central regions. 

The results of eq. (1191 and 1221 could be used for extending full 3-D simulations to densities up to n ~ 10 16 cm -3 . 
However, we have to remark that such an approach is likely to be correct only if the collapsing object remains approximately 
spherical at all stages, and if the spatial profiles (of density, temperature, chemical composition etc.) actually resemble the 
1-D results, at least until the moment when CIE cooling starts to overcome H 2 lines (n > 10 14 cm" 3 ). 



3 INSTABILITIES 

In this section, we investigate whether a collapsing protostellar cloud can become unstable to fragmentation by means of two 
different instability criteria. 
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Figure 6. Evolution of the ratio of the radiative cooling timescale t ra< j (thin solid line), chemical timescale £ c hem (dashed line) and 
cooling timescale t coo \ (thick solid line) to the dynamical timescale f d yn , as a function of density. The data come from the protostellar 
collapse simulations of R02 and refer to the central zone ("R02 path"). The three vertical dotted lines at densities ~ 3 X 10 13 , ~ 5 X 10 15 
and ~3x 10 16 mjj cm -3 show where the chemical heating term vanishes (implying t c hom = 00 an( i icool = *rad) an d changes sign. The 
horizontal dot-dashed line just visualizes the criterion t < t^ yn . 



3.1 Timescales comparison 

A classical criterion for fragmentation is that the cooling time t coo i must be shorter that the dynamical timescale 

• (23) 

In Fig. |S| we show the evolution of the ratios of three relevant timescales to the dynamical timescale, as determined from 
R02 data about the central evolution of a proto-stellar object (in the rest of this paper, we will refer to this evolutionary track 
as the "R02 path"). 

These timescales are: 

- the radiative cooling timescale (accounting for the effects of cooling) 

[l/(Mm H )]fc B r 

trad = 2W7-1) ( } 

where L ra d is the total radiative cooling rate per unit mass (L ra d = Lu ncs + I/cie in the conditions we are considering) and 7 
is the mean adiabatic index of the gas (cfr. ON98); 

- the thermo-chemical timescale (accounting for the thermal effects of chemical reactions, namely of formation or disruption 
of H 2 ) 

_ (n/ M )fc B T (n/^)k B T 
|-B chem |(7-1) |nH 2 |XH 2 (7-l) 

where E c hcm — — mh 2 Xh 2 is the rate of variation of chemical energy (per unit volume), hn 2 is the rate of H2 formation (per 
unit volume) and xh 2 = 4.48 eV is the binding energy of H2 ; 

- the "effective" cooling timescale (accounting for both chemical reactions and radiation): 

[l/( M m H )]fc B T 

Icool = : — -• (^D) 

(irad — -Bchcm/p)(7 ~ 1) 
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It can be seen in Fig.[S]that the ratio t coo \/tdyn is close to 1 for most of the considered evolution, although it never really 
crosses this critical value. Since the above criterion is only a sufficient condition for instability, we consider a more detailed 
instability criterion next. 



3.2 CIE induced instabilities 

3.2.1 Chemo -thermal instability criterion 

Silk (1983; see also Sabano & Yoshii 1977) considered a gas cloud in chemical and thermal equilibrium and investigates the 
conditions for the growth of isobaric perturbations. Assuming infinitesimal perturbations of the form T = Ti exp(ikx + cut), 
a dispersion relation 

Alu 2 + Blu + C = (27) 

is obtained, which can be used to assess the stability of a given gas cloud. Similar results were obtained by Saio & Yoshii 
(1986), and OY03, who were able to drop the assumption of chemical and thermal equilibrium. 

We use the values of the coefficients A, B and C which are given by equations (A19-A27) of OY03: they depend on the 
temperature T, on the density p, on the molecular mass fraction /h 2 , on the two functions T{p,T, /h 2 ) = ( rs ^ e °f change 
of the molecular mass fraction) and C(p, T, /h 2 ) (cooling rate per unit mass), and on their partial derivatives T p = dT/dp, 
T f = dT/dfH 2 , T T = dT/dT, C p = dC/dp, C f = dC/df H2 , and C T = dC/dT. 

Since A must be positive, it is possible to have growing perturbations (i.e. , Re(tj) > 0) if and only if C < 0, with (cfr. 
OY03, appendix A) 

C= -^(TC T -pC p -C)(T f + ^pT P +^ ; ) 

i 

6-/ H 



3=^(1 + ^tKTFt - P T P ) - ^(Tj + ± P T P + 3^ 

+ ^-i:j2 )t j TT T - P T p ), (28) 



where [i = (1 — /h 2 /2) 1 is the mean molecular weight (OY03 assumed a pure H-H 2 gas) and t« = [3tt/(32Gp)] 1/2 is the free 
fall timescale. 



3.2.2 The H2 formation and cooling rates (T and C) 

We now apply the above instability criterion (eq. I28|l to the density and temperature regime where optically thin CIE 
cooling could be dominant. For this reason, we will only consider the density range 10 13 cm -3 < n p < 10 17 cm" 3 (where 
n v — pX/m,H = nn + 2wh 2 is the number density of protons); we will also restrict ourselves to the temperature range 
400 K < T < 3000 K, because of the limitations of the approximate cooling rate we are going to assume (eq.|SJl. 

In this region of the density-temperature phase space, H2 formation and dissociation proceeds through the three-body 
reactions described by PSS83, so that the function T can be written as 

T = ^ = n/[2n*4(l - /h 2 ) 2 - k 5 fn 2 ] (29) 
with 

^ = n p (l-^), (30) 

while ki and k$ (the notation comes from PSS83) are the reaction rates for the formation and disruption of H2 , respectively. 
We use the PSS83 formation rate (fct), but we modify the dissociation rate (fcs) in order to better approximate the "high 
density" (n > 10 9 cm" 3 ) results of Martin, Schwarz & Mandy (1996): 

ki = 5.5xl0" 29 r 1 cm 6 s _1 (31) 

k 5 = 2.2xlO" 9 T a2 e"^(l-e~T) cm 3 s _1 (32) 

where B 5 = 51800 K, C 5 = 6000 K. 

At these high densities ON98 and R02 showed that we can safely assume that chemical equilibrium is attained, so that 
the condition T — can be coupled with equation 129^ in order to get the equilibrium H2 fraction fo(p,T): 

8n„fc 4 \ 1/2 " 



(33) 
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Figure 7. Value of C assuming no conversion of gravitational into thermal energy (<E> 



thin contours in the unstable region corresponding to C ■ 



-10" 12 , -10" 9 and -10" 6 s 



We show the n — T plane with the dashed 
the solid thick line for C = 0, and the thin 



dotted contours with C = 10 — 15 , 10 -12 and 10 -9 s — 2 in the stable part of the phase diagram. The thick dot-dashed line marks the "R02 
path" (the path followed by the central regions of the primordial proto— stellar cloud simulated in R02). The thin solid line, corresponding 
to /h 2 (Pi'^) = fo(p,T) = 0.5, separates the low temperature region where H is mostly molecular (and eo. Islcan be used) from the high 
temperature region where H is mostly atomic. 



However, we note that the third addendum in eq. 1)28^ introduces a strong dependence of the value of C upon the actual 
H2 abundance, and that even relatively small differences (a few percent) in /h 2 can lead to significant discrepancies in the 
results, as we will show in the next subsection. 

With regard to the "luminosity" £, in this subsection we will neglect all forms of radiative cooling except CIE optically 
thin cooling. In addition, the internal energy also changes because of the thermalization of gravitational energy. Following 
S83, we will assume that a constant fraction <E> of the gravitational energy in bulk motion is thermalized, so that 

£ ~ A C mpT a Xfu 2 - A*p 1/2 T (34) 



where 



.>/, 1 



A* = J^M=^ ~ 5.9 x 10 4 $ erg cm 1 ' 5 g" ' 5 KT 1 , (35) 



and we have expressed the CIE cooling rate using eq. Q. As we previously remarked, this last approximation requires the 
gas to be predominantly molecular (/h 2 ~ 0.5); in the regions of phase space where this is not true (that is, close to the upper 
temperature limit of 3000 K), the following results are probably inaccurate. 



3.2.3 The instability region 

We examine the sign of C across the n~T phase space by assuming that at each point we have an equilibrium H2 abundance 
(that is, fn 2 {p,T) — fo(p,T)) and that the values of T and L (and their derivatives) can be found through eqs. JU) and 
1341 1 . respectively (see the appendix for the explicit expressions of the derivatives). We still need to assume a value for the 
$ parameter in eq. dU, an d we choose to investigate two of the most relevant values: $ = and $ = 0.5. The former is 
consistent with the assumed spherical collapse scenario, while the latter is the value that can be expected if a disk in keplerian 
rotation forms, although such scenario is not consistent with the rest of our assumptions (for example, in a disk geometry the 
cooling rate is likely to be significantly higher). 

In Figures |7| and |H| we show contour plots of the values of C in the n — T plane in the two cases. 
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Figure 8. Same as Fig. 171 assuming a 50% conversion of gravitational into thermal energy ($ = 0.5). 



In both plots it can be seen that at high temperatures (T £ 1500 — 2500 K, depending on n and <E>) there is a region 
where C < and the chemo-thermal instability can occur. Since the conversion of gravitational energy into thermal energy 
has a stabilizing effect, the size of the unstable region keeps reducing with increasing values of but it never disappears, 
even when $ = 1. 



3.3 Instability on the R02 path 

In figures [7J and |H1 it is quite clear that the R02 path always remains outside the unstable region (<E> = 0.5), or barely grazes 
it ($ = 0). 

However, this results suffers from some important uncertainties: first of all, we already mentioned that the third addendum 
in eq. 1281 introduces a strong dependency of C on the exact value of the H2 abundance /h 2 (and of its time derivative J-); 
less importantly, it only considers optically thin continuum cooling, but this is not completely appropriate at the two extremes 
of the considered density range. 

Despite these subtleties the simple fact that the numerical results of R02 fall so close to the C = lines indicates 
that the presented analysis may be also applied to the collapsing proto-stellar cloud as a whole. Being able to describe the 
exact evolution and formation of the primordial proto-star purely analytically is highly desirable and will be the subject of a 
subsequent paper. 

A better investigation can be accomplished by applying the same criteria described in the previous subsections to the 
analysis of the behaviour of the object simulated by R02 (and also by ON98). Since we have access to the full evolution of 
such object, we can study its properties in greater detail, and in a larger range of densities (10 s cm -3 < n p < 10 17 cm -3 ). 

First of all, we can try to fix the problems about H2 abundance. This can be done in at least two different ways: 

(i) take /h 2 from R02 data (rather than by assuming /h 2 = fo) and leave everything else unchanged: in particular, we take 
the values of T and of its derivatives from eq. I|29ft ■ 

(ii) take /h 2 from R02 data, and also estimate values of T which are consistent with these data, but keep calculating the 
values of the partial derivatives Tt, T p and Tj from the analytical formulae we were using in the previous subsection. 

Approach (i) is the simplest, and it is very useful at low density, where the R02 abundances are very different from the 
equilibrium ones because the collapsing object had not reached chemical equilibrium yet, and where both ways of estimating 
T (from eq. 1291 and from R02 data) lead to identical results; on the other hand, at medium and high densities the results are 
very different (despite the small difference in the /h 2 values, generally just a few percent) , and it must be remarked that there 
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Figure 9. Value of the parameter C (top panels) and of the ratio of the fluctuation growth timescale 1/m to the free-fall timescale tg 
(bottom panels) along the R02 path, with two different assumptions about the H2 fraction variation rate T . In the left panels we take 
the H2 fraction /h 2 from the R02 results, but use the analytical values of T (case i). In the right panels we use the R02 data to get 
both /h 2 and tF, but we continue to estimate the derivatives of T from analytical formulae (case ii). In the upper panels, the ranges 
where C < (which implies instability) are denoted through a thick line, while thin dashed lines denote the ranges where C > (and 
the object is chemo-thermally stable). In the bottom panels we do not show this ratio when C > since that implies u> < 0. 



is a logical inconsistency in taking the H2 abundance and its rate of variation from two different sources. From this point 
of view, method (ii) is clearly superior, even if not optimal: from the purely logical point of view, we should obtain also the 
values of Tt, T p and T; from the R02 data, but unfortunately there is no way in which we can do that. 

Fixing the luminosity term is much simpler, since we don't have so many uncertainties: we add a term for the H2 line 
luminosity into the definition of the function C, and modify the CIE cooling rate by introducing a "CIE optical depth" 
correction factor. For the first task, we can employ the simple equation (1221 combined with the Hollenbach & McKee (1979) 
cooling rate, while for the second we have empirically found that a correction term of the form — - — - (with r c = (p/ po,c) 2 ' 8 and 
Po,c = 3.3 x 10 -8 gem -3 ), in conjunction with a slight reduction of the normalization (Acie,ro2 = 0.054 erg cm 3 g -2 s _1 K -4 , 
the difference is likely due to the different data set used by R02 for estimating CIE cooling) provides an excellent fit to the 
R02 CIE data up to the highest densities we are considering: the agreement with R02 data is actually good at least up to 
n ^ 10 19 cm -3 (figures 01 and |5] were obtained with this approximation). Finally, we will only consider the <& = case, so in 
the following we assume that 

C = ^Hmax[l, (p/pew)"' 3 '] + A C iE,no2T a Xpf ll2 1 ~ & - (36) 

The results are shown in Fig. El where we show the evolution of the value of the C parameter and of the corresponding 
instability growth timescale along the R02 path, as a function of density. The two sets of panels show the results obtained 
with each of the above assumptions about the values of /h 2 and T . 

First of all, it is apparent that the R02 object is unstable at low densities (n ,*$ 2 x 10 10 cm -3 ), and this is independent 
from the assumptions about T. This regime, where the instability is due to the fast H2 formation at the onset of 3-body 
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reactions, corresponds to the one found by PSS83 and S83, and recently discussed by ABN02 and OY03. However, the bottom 
panels clearly show that the instability growth timescale is always significantly longer than the free-fall timescale, so that 
this kind of instability cannot actually lead to fragmentation. This result is slightly different from the one recently obtained 
by OY03, who find that there exists a density range about n ~ 5 x 10 9 cm -3 where the growth timescale is slightly shorter 
than the free fall one, but their conclusion (the instability is present but does not lead to fragmentation) is the same as ours 3 . 
A different interpretation for the lack of fragmentation is the one pointed out by ABN02, that is, mixing due to turbulent 
motions: if these turbulent motions move at about the speed of sound c s , they should be able to erase any fluctuation of size 
^ c s /uj, and we find that the "core" where fluctuations could develop is always several times smaller than that. 

The second result is that CIE cooling can lead to an instability, at least if we estimate T from the R02 data (second 
approach): the instability occurs over about one order of magnitude in density, when continuum optically thin cooling is 
dominant. However, the bottom panels show that a situation similar to the one described at low density is very likely to 
occur: the instability growth timescale is always much longer than the free fall timescale, and fluctuations cannot be large 
enough to survive turbulent mixing. We also conducted some experiments taking reasonable guesses at the values of Tt, T p 
and J-f, with the result that the unstable range could extend over the whole phase when the cooling is due to optically thin 
continuum emission, but the fluctuation growth is always too slow to produce fragmentation. 



4 DISCUSSION AND CONCLUSIONS 

We have investigated the chemical and thermal instabilities of primordial metal-free gas during its collapse towards the 
formation of a protostar. We have focused on the Collision Induced Emission (CIE) and on its effects on fragmentation. 
We point out that even when H2 line cooling becomes optically thick at densities n Si 10 10 cm -3 this does not necessarily 
set a minimum mass scale for fragmentation due to an opacity limit. At higher densities (n Si 10 14 cm -3 ) the optically 
thin continuum (CIE) cooling dominates, and the collapsing protostellar cloud can once again fulfill the necessary (but not 
sufficient) fragmentation criterion t coo \ < idyn- We present accurate approximations to cooling and chemical rates that are 
particularily useful for analytical studies. A detailed analysis of the chemo-thermal stability of collapsing primordial gas 
leads us to conclude that the thermal instability arising from optically thin continuum cooling is of similar strength as the 
three body H2 formation instability discussed previously (PSS83; S83; OY03). However, our analytical analysis of the low 
density instability confirms the full 3D simulations results of Abel et al. (2002) in showing that formally the instability is 
very strong yet it cannot grow sufficiently fast to lead to independently collapsing fragments. In hindsight, it may not be too 
surprising that this instability does not lead to fragmentation. It fulfills the necessary (but not sufficient) condition of leading 
to a cooling time shorter than the dynamical time. An unstable patch of gas will, however, grow initially only iso-barically 
as it is compressed from the surrounding higher pressure regions. So although the cooling rate can increase dramatically as 
one increases the H2 fraction by a factor one thousand the minimum temperature the gas cools to cannot. Consequently, 
correspondingly small density fluctuations are formed which do not become gravitationally unstable and become mixed back 
in into the slightly warmer surrounding regions. 

Given the efficiency of collision induced emission cooling (Figure and its relatively small effect on the temperature in 
the one dimensional results one is tempted to conclude that it will also not lead to fragmentation. This is most likely a fair 
assessment since again the growth times are long and the initial iso-baric growth leads to only small overdensities. However, 
in a fully three-dimensional calculation the CIE cooling might be able to cool a disk forming around the primordial proto-star 
sufficiently far in order to become gravitationally unstable and form stellar, or even planetary size companions (see e.g. Boss 
1993 and Boss 2002). Given the dramatic accretion rates of these proto-stars, however, it is not clear whether a sufficiently 
massive disk may form in the first place preventing a conclusive answer. 

Despite the simplicity of our derivation the analytical rates for the H2 line cooling rate in the optically thick regime 
and optically thick continuum cooling they are in remarkably good agreement with the radiative transfer calculations of R02 
and ON98 that followed the detailed transfer of hundreds of molecular hydrogen rotational and vibrational lines. We have 
shown that implementing these simple approximations into the one dimensional Lagrangian hydrodynamics code developed 
by R02 leads to virtually indistinguishable results over all relevant density and temperature ranges. This clearly demonstrates 
that at least in the initial phases of primordial proto-stellar formation the effects of radiative transfer are in essence a local 
correction to cooling rather than a transport of energy to distant regions of the hydrodynamic flow. This is an important result 
which, given the analytical and purely local cooling correction factors derived here should one allow to follow the proto-stellar 
collapse in three dimensions all the way to stellar densities employing the recently developed versions of 128bit adaptive mesh 
refinement techniques of Bryan, Abel & Norman (2001) or similar high dynamic range hydrodynamic codes. 



3 This difference couid partiaiiy arise from the slightly different definition of dynamical timescale used by OY03; see their eqn. A16. 
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APPENDIX A: SUMMARY OF FORMULAE AND DERIVATIVES 
Al Luminosity 

The luminosity per unit mass C can be written as (see eq. 1361 

£ ~ ^lHmax[l, (p/po.i)"' 3 '] + A cm ,nmT a Xpfn 2 — — - A^Tp 1/2 , (Al) 

where we have included the thermalization of gravitational energy term. The constants values are po,i = 1-34 x 10~ 14 gem -3 , 
Pi = 0.45, A C ie, R02 = 0.054 ergcm 3 g- 2 K- 4 , a = 4, r c = (p/p , c ) 2 ' 8 (po,c = 3.3 x 10~ 8 gem" 3 ) and A* = 5.9 x 
10 4 erg cm 1 ' 5 g" ' 5 K" 1 ; H(T 3 ) is given by eq. 0. 

If the collapsing object remains approximately spherical, this formula is reasonably accurate up to densities p ~ 2 x 
10 -5 gem" 3 . 

Because of the maximum in the H2 lines term, it is useful to distinguish two cases, 
if p < po,i, the partial derivatives of C are 



Ct = f§= -7^^f +aAcmT X P-f^ ~ A *P ! (A2) 



aAcml Xpfn 2 A$p ' 

1 _ o-T'c 1 1 

o£ — AciET a Xfn 



m-H 01 t, 
1.... v-vf.. 3i . e + _ 3<; )JjliL 



\A*Tp- x l 2 (A3) 
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Cf = 



X 
m H 



1 -e" 
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(A4) 



Instead, if po,i < p we have: 



x fa 2 an 

m H ST 



*-H(p/po,i)- pi + A C mT a Xp^ 



where, in both cases, 



(A5) 
(A6) 
(A7) 



dH dT 3 on icr 



dT dT dT s T\ 



9.5 x 10~ 22 r| 76 ^ e -(%P) 3 f 3.76T 3 +0.2T 3 3 1 + 0.0066 ^ 



1 + 0.12T, 



2.1 

3 



1 + 0.12T. 



2.1 



T 2 

3 



+ 



10" 



1.5 x 10" 2 e T 3 + 3.9 x 10" 18 e T 3 + 1.9 x 10" 17 e T s 



ergs K 



(A8) 



A2 H 2 Formation 

The H2 formation function is given by 



dt 



n/[2npfe 4 (l - /h 2 ) -fc 5 /H 2 ] 



where 

n p = pX/m H , nf = np(l - 15/h 2 /16) 
and (as explained in section 3.2.2) 

fc 4 = A4T" 1 , k 5 = A 5 r°- 2 e - S5/T (i - e -° 5/T ), 

with A A = 5.5 x 10~ 29 cm 6 s" 1 , A = 2.2 x 10" 9 cm 3 s" 1 , B 5 = 51800 K, C 5 
The partial derivatives are 

1 



6000 K. 



dT 
dT 



nj_ 

T 



r= %= ^[4n p fc4(l-/H 2 ) 2 -fc 5 /H 2 ] 



2n p fc 4 (l — /h 2 ) 2 + k 5 fn 2 ( 

2 



no ^ 

T T e^/T-i 



15 



16 - 15/ 



[2n p fc 4 (l - /h 2 ) 2 - fc 5 /H 2 ] + 4n p fc 4 (l - /h 2 ) + h 



(A9) 
(A10) 
(All) 

(A12) 
(A13) 

(A14) 
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